%% Compare Plain VARs and Panel VARs
% by Jaromir Benes
% 
% Compare a number of reduced-form and structural properties of the VARs
% estimated on individual countries, and the VAR estimated on a panel of
% data. Show the eigenvalues, the implied unconditional cross-correlation
% coefficients, estimated residuals, and impulse responses from structural
% VARs (using a simple Cholesky identification scheme).

%% Clear Workspace
%
% Clear workspace, close all figure windows, move to the top of the command
% screen, and check the IRIS version.

clear;
close all
home;
irisrequired 20130401;

%% Load Individual Plain VARs and Panel VAR
%
% Load the VARs estimated on individual countries, `vau`, `vca`, and `vno`,
% and their output data, `dvi`, from the mat file created in
% `estimate_individual_countries`. Load the estimated panel VAR, `vp`, and
% its output data, `dvp` from the mat file created in `estimate_panel`. You
% must run `estimate_individual_countries` and `estimate_panel` at least
% once before this m-file to create these mat files.

load estimate_individual_countries.mat vau vca vno dvi;
load estimate_panel_var.mat vp dvp;

%% Calculate Unconditional Covariances/Correlations
%
% Calculate the implied unconditional cross-covariances and
% cross-correlations up to first order, both from the VARs estimated on
% individual countries, and the panel VAR. All the cross-covariance and
% cross-correlation matrices have the same size <?size?>: Ny-by-Ny-by-K,
% where Ny is the number of variables included in the VAR, and K is the
% order up to which the covariances and correlations are computed (NB: this
% is not the order of the VAR!).

[Cau,Rau] = acf(vau,'order=',1);
[Cca,Rca] = acf(vca,'order=',1);
[Cno,Rno] = acf(vno,'order=',1);

[Cp,Rp] = acf(vp,'order=',1);

size(Cau) %?size?

%% Compare Cross-Correlations
%
% Use the function `plotmat` to visualise the contemporaneous
% cross-correlatins. The function depicts the matrix entries as disks with
% their diameters proportional to the size of the entry (i.e. the
% cross-correlation). Positive entries are blue, negative entries are red.
% Use the option `'showDiag=' false` <?showdiag?> to hide the diagonal
% entries; they always equal `1` in contemporaneous cross-correlation
% matrices, and there is no need to display them.

figure();

subplot(2,2,1);
plotmat(Rp(:,:,1),'showDiag=',false); %?showdiag?
grid on;
title('Panel');

subplot(2,2,2);
plotmat(Rau(:,:,1),'showDiag=',false);
grid on;
title('Australia');

subplot(2,2,3);
plotmat(Rca(:,:,1),'showDiag=',false);
grid on;
title('Canada');

subplot(2,2,4);
plotmat(Rno(:,:,1),'showDiag=',false);
grid on;
title('Norway');

grfun.ftitle('Cross-correlation coefficients');

grfun.ftitle('off-diagonal elements only, blue=positive, red=negative', ...
    'location=','south');

%% Compare First-Order Autocorrelations
%
% The implied unconditional first-order autocorrelations are found on the
% diagonal of the first-order cross-correlation matrix, <?diag?>. Plot the
% autocorrelations from the VARs estimated on individual countries (blue)
% against the those implied by the panel VAR (red) as bar graphs.
% Obviously, the red bars are the same in each of the graphs.

rhoau = diag(Rau(:,:,2)); %?diag?
rhoca = diag(Rca(:,:,2));
rhono = diag(Rno(:,:,2));
rhop = diag(Rp(:,:,2));

yList = get(vp,'yList');

figure();

subplot(2,2,1);
bar([rhoau,rhop]);
grid on;
axis tight;
set(gca(),'xTickLabel',yList,'yLim',[0,1]);
title('Australia');

subplot(2,2,2);
bar([rhoca,rhop]);
grid on;
axis tight;
set(gca(),'xTickLabel',yList,'yLim',[0,1]);
title('Canada');

subplot(2,2,3);
bar([rhono,rhop]);
grid on;
axis tight;
set(gca(),'xTickLabel',yList,'yLim',[0,1]);
title('Norway');

grfun.ftitle('First-order autocorrelation');

grfun.bottomlegend('Country VAR','Panel VAR');

%% Compare Eigenvalues
%
% Get the eigenvalues from the VAR objects estimated on individual
% countries, and those from the panel VAR, and plot them all in one graph
% against the unit circle. The eigenvalues determine the properties of
% structural VAR impulse (shock) responses.

eau = eig(vau);
eca = eig(vca);
eno = eig(vno);
ep = eig(vp);

figure();
hold all;
ploteig(eau,'unitCircle=',false,'quadrants=',false,'marker=','v');
ploteig(eca,'unitCircle=',false,'quadrants=',false,'marker=','o');
ploteig(eno,'unitCircle=',false,'quadrants=',false,'marker=','x');
ploteig(ep,'marker=','*');
axis tight;
grid on;
legend('Australia','Canada','Norway','Panel');
title('Eigenvalues');

%% Compare Estimated Residuals

eList = get(vp,'eList');
countryList = get(vp,'groupNames');

for iCountry = 1 : 3
    figure();
    country = countryList{iCountry};
    for iName = 1 : 4
        name = eList{iName};
        subplot(2,2,iName);
        h = plot([dvi.(country).(name),dvp.(country).(name)]);
        axis tight;
        grid on;
        title(['Residuals ',name],'interpreter','none');
    end
    grfun.ftitle(['Estimated residuals ',country]);
    grfun.bottomlegend('Individual VAR','Panel');
end

%% Compare Shock Responses
%
% Compare shock (impulse) responses from the VARs for individual countries
% and the panel VAR. Shock responses can only be calculated from
% structural VARs (SVARs).
%
% * <?svar?> First, convert the estimated reduced-form VARS to structural
% VARs (SVARs) by using the default Cholesky identification scheme, i.e. a
% lower triangular matrix of instantaneous responses. See `help SVAR/SVAR`
% for other identification schemes and options available.
%
% * <?srf?> Run the shock response function, `srf`, to get databases with
% the responses of all variables to all shocks.
%
% * <?graph?> Cycle over the shock response databases, and plot each
% response for all countries and the panel in one graph. Use the list of
% VAR variables, `yList`, <?yList?> to label the graphs. Each entry in the
% database is a multivariate tseries object with individual columns being
% the response in the respective variable to individual shocks <?columns?>.

svp = SVAR(vp) %#ok<NOPTS> %?svar?
svau = SVAR(vau) %#ok<NOPTS>
svca = SVAR(vca) %#ok<NOPTS>
svno = SVAR(vno) %#ok<NOPTS>

sp = srf(svp,1:20); %?srf?
sau = srf(svau,1:20);
sca = srf(svca,1:20);
sno = srf(svno,1:20);

yList = get(vp,'yList'); %?yList?

for iShock = 1 : 4 %?graph?
    figure();
    for iName = 1 : 4
        name = yList{iName};
        subplot(2,2,iName);
        h = plot([ ...
            sp.(name){:,iShock}, ... %?colums?
            sau.(name){:,iShock}, ...
            sca.(name){:,iShock}, ...
            sno.(name){:,iShock}, ...
            ]);
        set(h(1),'linewidth',2);
        axis tight;
        grid on;
        title(name);
    end
    grfun.ftitle( ...
        sprintf('Response to structural shock #%g',iShock), ...
        'interpreter','none');
    grfun.bottomlegend('Panel','Australia','Canada','Norway');
end

%% Help on IRIS Functions Used in This M-file
%
% Use either `help` to display help in the command window, or `idoc`
% to display help in an HTML browser window.
%
%    help VAR/Contents
%    help VAR/acf
%    help VAR/eig
%    help VAR/get
%    help SVAR/Contents
%    help SVAR/SVAR
%    help SVAR/srf
%    help grfun/bottomlegend
%    help grfun/ploteig
%    help grfun/plotmat
%    help grfun/ftitle